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■ ABSTRACT 

<D ■ 

Context. Most of our knowledge of the nonlinear development of the magnetorotational instability (MRI) relies on the 
QQ . results of numerical simulations employing the shearing box (SB) approximation. A number of difficulties arising from 

' this approach have recently been pointed out in the literature. 

' Aims. We thoroughly examine the effects of the assumptions made and numerical techniques employed in SB simula- 

I— I , tions. This is done to clarify and gain a better understanding of those difficulties, in addition to a number of additional 

■ serious problems raised here for the first time, and of their impact on the results. 

$ZL(' Methods. We used analytical derivations and estimates as well as comparative analysis to methods used in the numerical 

study of turbulence. We performed numerical experiments to support some of our claims and conjectures. 
Results. The following problems, arising from the (virtually exclusive) use of SB simulations as a tool for the under- 
standing and quantification of the nonlinear MRI development in disks, are analyzed and discussed: (i) inconsistencies 
^ ' in the application of the SB approximation itself; (ii) the limited spatial scale of the SB; (iii) the lack of convergence of 

1 I ' most ideal mgnetohydrodynamical (MHD) simulations; (iv) side-effects of the SB symmetry and the non-trivial nature 

of the linear MRI; and (v) physical artifacts arising from the very small box scale due to periodic boundary conditions. 
. Conclusions. The computational and theoretical challenge posed by the MHD turbulence problem in accretion disks 

^ ' cannot be met by the SB approximation, as it has been used to date. A new strategy to confront this challenge is 

' proposed, based on techniques widely used in numerical studies of turbulent flows - developing (e.g., with the help of 

I local numerical studies) a sub-grid turbulence model and implementing it in global calculations. 
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1^ 1. Introduction ployed to determine its value. Such an approach also in- 

' , , r volves identifying an instability (or combination of insta- 

t-^ ■ A physical mechanism for angular naomentum transport ^-^^-^^^ ^j^^^ responsible for the emergence of accretion 



>< 



outward is needed for the theoretical modelling of accre- ^-^^ turbulence. This quest, according to this strategy, has 
tion disks. Prendergast & Burbidge (1968) proposed the ^^^^^^ ^^^^^ -^^^^^^ ^^^^ sometimes contentious) 



^ existence of thm accretion disks to explain some galactic research for the three decades following the above men^ 

^ . X-ray sources. They realized that such disks must be char- ^-^^^^ ground-breaking studies. Notwithstanding the cur- 

acterized by an enhanced viscosity, orders of magnitude ^^^^ ^^^j^ ^ g^^^^^l ^^^^^^ turbulence, the identifi- 

larger that the microscopic viscosity expected for the gas ^^^j^^ instability process and, hopefully, also un- 

comprising the disks. Given the likelihood that such flows jerstanding the transition from laminar flow into some 

are turbulent, they estimated the value of this enhanced turbulent-like steady state (or at least quasi-steady one) 

viscosity with the help of a mixmg-length theory. A few ^^^^^^ ^^^^ promising strategy to quan- 

years later Shakura & Sunayev (1973) and Lynden-Bell & ^^g^j^^ momentum transport and energy generation 

Prmgle 1974) successfully bypassed the specific lack of un- accretion disks. However, no definitive candidate hydro- 

derstandmg of turbulence (a situation that is largely still dynamical instability had been convincingly identified and 

with us) by employing physically motivated parametriza- ^-^^^^ accepted before Balbus & Hawley (1991) discovered 

tions of this turbulent viscosity In doing so thm accretion ^^^^ ^^^^ magnetic fields destabilize differential Keplerian 

disks, for the first time, could be modeled m a variety of .station. The hnear magnetorotational instabihty (MRI; 

astrophysical settings. In turn, this has brought about sig- Velikhov 1959, Chandrasekhar 1960) was shown to oper- 

nificant progress m our understanding of these important ^^^^^ conditions characterizing rotationally-supported 

systems. ^ magnetized cylindrical disks. 

Because the value and nature of the now famous nondi- 

mensional viscosity parameter a still remains a theoreti- This important study has inspired intensive research ac- 

cal unknown, phenomenological approaches are usually em- ti^ity on the^ linear MRI and its nonlinear development. 

This mechanism has been investigated for various physi- 
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matical methods. The complexity of its nonlinear develop- 
ment requires employing numerical computational meth- 
ods. Understanding, or at least numerically simulating in 
a faithful way, the nonlinear transition and saturation pro- 
cesses (even through the use of some approximations), is 
indispensable if the aim is to improve the more than 30 
years old phenomenological a-model. 

As the general problem of the development of magne- 
tohydrodynmical (MHD) turbulence in accretion disks is a 
formidable one, various approximations have been utilized. 
The purpose of this paper is to assess the viability of one 
of the approximations most commonly used in numerical 
studies - the shearing box (SB), also known as the "shear- 
ing sheet." The review articles by Balbus & Hawley (1998) 
and Balbus (2003) contain a comprehensive summary of 
the results of linear studies as well as SB numerical ones 
(under various physical conditions) available at the time of 
these reviews, together with a thorough list of references, 
of which only some will be referred to here. 

In the recent years, a number of new and improved (in 
resolution) SB calculations subject to differing conditions 
and employing different numerical methods have been re- 
ported on. Some results of these studies have raised a num- 
ber of significant issues regarding a central question per- 
taining to accretion disks: a quantitative estimate and scal- 
ing of the angular transport in a saturated MHD turbulent 
state (c.g, Pessah, Chan & Psaltis 2007; Lesur & Longaretti 
2007, Fromang et al. 2007). It is of no surprise that the es- 
timate of this transport, or the "effective a", which is a 
quantitative measure thereof, must depend upon those rel- 
evant nondimensional numbers characterizing the flow. In 
particular, these numerical studies indicate that the trans- 
port is proportional to some positive power of the magnetic 
Prandtl number of the flow. 

In parallel, a number of related approximate analytical 
treatments using asymptotic analysis have been published 
(Knobloch & Julien 2005; Umurhan, Menou & Regev 2007; 
Umurhan, Regev & Menou 2007, hereafter URM07). The 
latter two studies found the above-mentioned scalings in 
some limits of a magnetic Taylor- Couctte thin gap config- 
uration (formally different from the SB only in the bound- 
ary conditions), while the former work showed that trans- 
port diminishes in proportion to the inverse of the Reynolds 
number and/or the magnetic Reynolds number for the non- 
linear development of coupled channel mode (CM) solu- 
tions. 

Recent years have also seen experimental efforts to ob- 
serve the MRI and its nonlinear development in the lab- 
oratory. These have been done primarily in a magnetic 
Taylor-Couette (mTC) setup with an applied axial or heli- 
cal magnetic field. In these studies (see, e.g., Ji, Goodman 
& Kageyama 2001; Noguchi et al. 2002; Sisan et al. 2004; 
HoUerbach & Riidiger 2005; Liu, Goodman & .li 2006; 
Stefani et al. 2006; and references therein) MRI turbulence 
has only been reported for experiments with helical fields, 
in which the azimuthal field component is significant. Even 
for this case some controversy remains over the interpreta- 
tion of the results (Liu et al 2006; Riidiger & HoUerbach 
2007). 

Most of the existing SB calculations (at least of the 
early ones) employed a base state with a constant verti- 
cal magnetic field and followed the development of pertur- 
bations on this state. This situation, usually referred to 
as non-zero, or fixed net flux (herafter fixed-flux), seems 



to be the most appropriate for disks threaded by exter- 
nal magnetic flelds (whether primordial or originating from 
the central accreting object). The majority of existing nu- 
merical works of this sort have used ideal (that is, non- 
dissipative) MHD SB equations with periodic boundary con- 
ditions (periodic-BC) in the manner employed, for example, 
by Hawley et al. (1995), but the very recent calculations of 
Lesur & Longaretti (2007) include explicit viscosity v and 
resistivity 77 in a small shearing box (SSB, see below), high- 
resolution, numerical simulation. Systematically changing 
the Reynolds number (Re) and the magnetic Prandtl num- 
ber Pm = I'/f], in some range, these authors aimed at un- 
covering the trends and scalings of relevant physical prop- 
erties, notably angular momentum transport, with the rele- 
vant non-dimensional number (s). A number of studies also 
exist demonstrating the role of the MRI as an essential in- 
gredient for internal dynamo action in the case where the 
base state has zero net flux (hereafter zero-flux) . It is rea- 
sonable to suggest that the zero-flux case shoiild be relevant 
to disks (or disk regions) that are not threaded by magnetic 
flux. The most recent zero-flux calculations in the SB frame- 
work with relatively high resolution have been reported by 
Fromang & Papaloizou (2007) and by Fromang et al. (2007) 
These studies (both employing periodic-BC) focus on the 
issue of convergence (i.e., numerical resolution) in ideal SB 
calculations and on the Re and Pm dependencies of the 
transport in dissipative conditions. 

Recently King, Pringle & Livio (2007), hereafter 
KPL07, pointed out the discrepancy (of at least one order of 
magnitude) between the values of a inferred from observa- 
tions and the estimates of this parameter based on (mainly 
the early) SB numerical simulations. In particular they re- 
examined the conclusion, reasoned by Hawley, Gammie & 
Balbus (1995) based on their flxed-flux SB simulations, that 
the effective a is dependent (almost linearly) on the value 
of the imposed Bz- Similar scaling has also been reported 
by Sano et al. 2004 and Pessah, Chan & Psaltis 2007, but 
the latter also found a pre-factor depending on the box 
size as L^-^. They concluded from this puzzling result that 
fixed-flux simulations are not likely to be adequate for ac- 
cretion disks and the zero-flux ones should be used instead. 
They also identified several theoretical difficulties of the 
SB simulations - limitation of scale, problematic boundary 
conditions, question of convergence of the simulations (nu- 
merical resolution issues) , and a possible breakdown of the 
MHD approximation. Rincon, Ogilvie & Proctor (2007) re- 
cently published a significant study, hereafter ROP07, in 
which they discovered that a self-sustaining nonlinear dy- 
namo process may operate in Keplerian shear flows in the 
zero-flux case. In this dynamo process, the MRI is not a 
bona-fide linear instability, resulting from a modal stabil- 
ity analysis of the base flow, but rather participates in one 
of the stages of the dynamo action. We shall refer to the 
role of the MRI in this process as mediating rather than 
"driving" , or "inducing" (which will be the term in the 
fixed-flux case) MHD turbulence. It is important to note, 
in the context of the present work, that this analysis was 
done in the framework of a rotating magnetic plane Couette 
flow (rmpC), that is, one in which realistic wall boundary 
conditions (wall-BC) were employed. This is in contrast to 
the majority of SB calculations wherein periodic-BC are 
used instead. 

Tiirbulent magnetic dynamo theory, which is concerned 
with the question of magnetic fleld generation in cosmic 
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bodies (see, e.g., Moffat 1978; Parker 1979), is a principally 
similar MHD problem and it has been studied extensively 
using direct numerical simulations (e.g, Cattaneo & Hughes 
1996; Brandenburg & Donner 1997; Schekochihin et al. 2007 
and references therein). The issues regarding the computa- 
tional domain, resolution and BC in such calculations, and 
their effect on the results and findings have always seri- 
ously been considered (see, e.g., Blackman & Field 2000). 
We mention this theory here at the outset because we will 
use it as an example of another intensively studied MHD 
turbulence problem in which progress has been quite im- 
pressive. 

In this paper, we shall examine in detail the viability of 
the SB approximation to the study of angular momentum 
transport (resulting from MHD turbulence related to the 
MRI) in accretion disks. We shall discuss it for both the 
fixed-flux case (where the MRI is expected to induce MHD 
turbulence via a supercritical transition, i.e., it is a linear 
instability in the usual modal sense) and the zero-flux case 
(where the MRI is invoked in a self-sustaining process in- 
volving a subcritical transition) . We shall not only address 
and expound on some of the issues discussed by KPL07 and 
others, but also point out and analyze in detail a number 
of other very significant difhculties. This will lead us to the 
conclusion that the SB approximation, as used in the great 
majority of the existing numerical simulations, suffers from 
limitations that make it inadequate for the qualitative and 
quantitative understanding of the nonlinear development of 
the MRI in accretion disks. 

Our paper is organized as follows. We start by review- 
ing a recent fairly exact derivation of the SB approximation 
distinguishing between the two limits of this approxima- 
tion: the small shearing box (SSB) and the large shearing 
box (LSB). In the next section the relevant physical length 
scales of the problem are discussed and related to the box 
scale, yielding an answer to the question: which physical 
properties can be expected to be faithfully uncovered in 
SSB and in LSB? In Sect. 3] we consider numerical resolu- 
tion (convergence) in SB simulations and its implications on 
the numerical results and the conclusions thereof. In Sect.[5l 
we point out the symmetry properties of the SB equations 
and of the boundary conditions used in their simulations 
(which are not globally valid in an accretion disk model) 
and their effect on the solutions to the linear and nonlinear 
problems. This is followed in Sect. [H] by a discussion of the 
SB energetics and the effects of periodic boundary condi- 
tions on it, when the box size is not large enough. We also 
perform a number of numerical experiments to demonstrate 
this important issue. 

As far as investigating the properties of MRI induced 
and/or mediated MHD turbulence is concerned, wc finally 
discuss which (if any) of the setups - SB, mTC or rmpC 
- would be an appropriate venue for local numerical stud- 
ies. In this context, we stress the value of analytical and 
semi-analytical studies for the purpose of better under- 
standing numerical simulations and experiments. We also 
suggest possible promising strategies for conducting global 
(disk scale) studies, exploiting the knowledge gained from 
the local ones. The ultimate goal in all of this is to find a 
viable way to quantitatively assess the turbulent angular 
momentum transport and energy dissipation for the mod- 
eling of accretion disks. 



2. The Shearing Box (SB) approximation 

The essence of the SB approximation is in its local approach, 
that is, the resulting equations for the perturbations on 
a steady base flow are approximately valid in a small re- 
gion (a Cartesian box) around a typical point in the disk. 
The global MHD base flow in an almost Keplerian accre- 
tion disk is not only rotating and strongly sheared, but it is 
also inhomogeneous, non- isotropic (endowed with non-zero 
gradients in density and other physical variables and these 
gradients have very different scales in different directions) 
and swirling (streamlines are curved). 

The SB approximation, with its emphasis on the effect 
of shear, is not new and has been introduced in an as- 
trophysical context long ago by Goldreich & Lynden Bell 
(1965) (see also Toomre 1964) in their study of sheared 
gravitational instabilities in a galaxy. Hawley & Balbus 
(1991) adapted the approximation to the numerical study 
of the MRI in accretion disks and it has been routinely used 
ever since. Umurhan & Regev (2004) followed a systematic 
derivation of the approximation, using asymptotic scaling 
arguments with the purpose of quantifying the approxi- 
mations that are made leading to the SB (see Appendix 
A of that paper). Even though a purely hydrodynamical 
case was considered there, the addition of MHD terms is 
straightforward, because the magnetic field remains invari- 
ant in a transformation to a moving frame as long as the 
velocities involved are non-relativistic (naturally with re- 
spect to an observer in the laboratory frame) . 

We summarize here the procedure for an MHD base flow 
in a thin Keplerian disk with a constant vertical magnetic 
field (as a convenient example). This is done to set the 
stage and also emphasize the difference between the two 
self-consistent limits of the approximation. The derivation 
starts by picking a point in the disk, (tq, ^ = 0) in cylin- 
drical coordinates, transforming the relevant full equation 
set into a frame rotating with the Keplerian angular veloc- 
ity appropriate for this point = ^^k(''o)7 f^nd looking at 
the equations in a Cartesian box - a small neighborhood 
around the point tq. Non-dimensionalization of the equa- 
tions that describe a thin disk brings out two small non- 
dimensional parameters - e (the typical disk height /ig, in 
units of ro) and S (the box size in the same units). The next 
step is the expansion (up to first order in 6) of the depen- 
dent variables of the laminar (mean) base flow, Keplerian 
rotation in this case, in the neighborhood of tq, and the re- 
moval of the mean flow. This results in a set of partial differ- 
ential equations for the disturbances that are fluctuations 
atop a base state (i.e, deviations from a laminar flow pro- 
file with constant vertical magnetic field). Identifying the 
radial (shear-wise) direction with the Cartesian coordinate 
X, the azimuthal (stream- wise) direction with y, leaving z as 
the vertical coordinate and neglecting the curvature terms 
(consistently with the smallness of 6) finally yields the SB 
equations. For details see Umurhan & Regev (2004), who 
identified two limits of the SB approximation, depending 
on the ordering of the small parameters e and S. In the first 
one, a box whose size is much smaller than the vertical scale 
height of the disk, i.e., 6 e is considered, and thus the 
unperturbed state of linear shear may be considered ho- 
mogeneous. The perturbations are then incompressible and 
acoustics are filtered out. This is the SSB. The second limit, 
LSB, is the case i5 = O (e), in which vertical stratification 
as well as compressibility must be taken into account and 
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the situation is considerably more comphcated than SSB. 
We note here in passing that a typical thin accretion disk 
has e ~ 10~^, and this is thus the LSB {S = e) scale in units 
of tq, while the SSB {6 = e^) scale is in these units ~ 10^*. 

For the purposes of this paper, it will suffice to con- 
sider here the simple SSB for dissipative MHD flow with a 
fixed vertical background field, possibly zero, and constant 
density (= 1 in our units), because the issues we want to 
deal with do not depend on additional details. The velocity 
perturbations (deviations from the laminar profile, in units 
of floToS) are u = {u,v,w) and those in the magnetic field 
(expressed in units of V^, a typical Alfven speed squared, 
appropriate for the background field Bq) are correspond- 
ingly b = (bx,by,bz). The total pressure perturbation (in 
units of e^Slp^'oi see below) is written as zu = p + (2/3)~^&^, 
where /3 = e^ri^rg/yj (the "plasma beta"). Additionally, 
lengths are scaled by the box size (Sro) and time by the 
Keplerian time (l/Jlo). With these assumptions and defi- 
nitions the non-dimensional SSB equations in the rotating 
Cartesian frame explained above read, 

V-u = 0, V-b = 0, (1) 



Keplerian velocity is highly supersonic. As we have pointed 
out above, on the SSB scale the fiow is essentially incom- 
pressible and the sound speed loses its meaning. By scaling 
the pressure with e^rigrp, we thus consistently express the 
vertical size of the disk. 

It is possible to reformulate the above SSB equations 
in terms of shearing coordinates in the same way as imple- 
mented by Goldreich & Lynden-Bell (1965). These coordi- 
nates are written into a frame that is shearing exactly as 
the background fiow itself on the SB scale (linearly). This 
coordinate system (^, 77, (, r) is thus defined in terms of the 
coordinates in the rotating frame through the transforma- 
tion 

£,^x, ri = y-qnotx, ( = z, t = t, (6) 

and then equations ([TJS]) may be transformed in terms of 
the shearing coordinates by only making the following for- 
mal replacements 

dt dr-q^Q^d,,; dx d^^-qVloTdn] dy 9,,, dc_{l) 
and 



(dt — qQ.Qxdy)u — 2Q.qv + u • Vu = —dx'ou + 

+ ^ (Bo a, + b • V) 6, + ^V^u, (2) 
p Ke 

(dt — qfloxdy)v + (2 — q)^ou + u • Vw = —dyW + 

+ \{Bod,+h-V)hy + ^V^w, (3) 
p Ke 

{dt — qfloxdy)w + u • Vw = —dzzu + 

+ ^ (Bo a, + b • V) fe, + -^y^w, (4) 

p Ke 

{dt - qil.oxdy)h = BodzU + 

+ qilob^y + V X (u X b) + -^V^b, (5) 

where the Reynolds number (Re) and its magnetic coun- 
terpart, Rcm = Re • Pm, are defined in the SSB as: 

Re.Mf, Re„.^^. 
V rj 

For the sake of some generality and convenience we have 
kept r^o and Bq (which are actually equal to 1 in our units) 
and q = — ((ilnri/(ilnr)o (equal to 3/2 for a Keplerian ro- 
tational law). Note that to the extent to which these equa- 
tions represent a small disk section, if the base fiow devi- 
ates significantly from Keplerian rotation, then the hori- 
zontal pressure gradient does not drop out of the equations 
and the system would not be a proper representation of a 
rotationally-supported disk. In addition, it is perhaps im- 
portant to explain the choice of the pressure unit. In a disk 
that is vertically supported by pressure, the vertical scale- 
height is determined by the sound speed and gravity (that 
can be expressed in terms of the local Keplerian speed). 
Thus at a point ro , we obtain that this vertical scale is ap- 
proximately given by — Vs/^lQ, where Vs is the typical 
sound speed. Therefore e = Vs/ {^^fa) find it is small if the 



V->0-^gr!or5^, (8) 

where D = -I- f}dri + C^Q is the gradient operator in the 
shearing coordinate system. 

The LSB equations (both in terms of an observer in 
the rotating frame and in terms of the shearing coordi- 
nates) are somewhat more complicated because the base 
state (background) vertical profiles explicitly appear and 
the fiuid must be treated as compressible. Thus, equations 
of mass as well as energy conservation have to be included, 
together with some kind of constitutive thermodynamic clo- 
sure relation (see Umurhan & Regev 2004, for details). If 
the LSB approximation is to be used in the context of an 
accretion disk (this kind of approach is sometimes called 
quasi- global), the vertical BC are no longer a matter of 
"free" choice as they are actually the physical conditions 
on the vertical edges of the disk. 

We remark here, at the outset, that a significant number 
of published SB numerical simulations seem to use neither 
SSB nor LSB, but a kind of "intermediate" approximation 
between them. As the SSB and LSB are asymptotic repre- 
sentations of thin-disks, it is not clear how much and which 
of the results derived from these intermediate approxima- 
tions are physically relevant to disk systems. A number 
of examples of such seemingly inconsistent calculations are 
given in the next section. 

We finally note that virtually all numerical calculations 
that utilize the SB approximation, in all its variants, employ 
periodic-BC in which the periodicity in the radial direction 
is sheared (see Hawley et al. 1995). These sheared-periodic- 
BC are equivalent to enforcing that all perturbation quan- 
tities be triply (or in the case of a 2-D calculation, doubly) 
periodic in the sheared coordinates system. We shall discuss 
all these issues and their implications in more detail below. 

3. Relevant length scales 

In a typical accretion disk, ro(= 1 in units of ro, which we 
shall use in this discussion as the length unit) is the scale 
on which curvature terms appear and this is also the scale 
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of any underlying radial structure gradients. Vertical strat- 
ification obviously appears on the vertical pressure scale- 
lieight, which is e = Hq/tq <C 1 in our notation and units. 
This is also the scale on which the effects of compressibility 
are manifested. These are the basic physical disk structure 
length scales. 

In studies of turbulence it is customary to consider the 
injection scale, i = , on which the process of energy in- 
jection into the system is effected and the dissipation scale, 
= k'^^ , on which microscopic viscosity (and in MHD also 
resistivity) are operative. These are linked by the inertial 
range of scales. Since fcin <C fed (at least in astrophysical 
systems), the inertial range extends over many orders of 
magnitude. In MHD turbulence, the energy cascades in 2D 
and 3D are both direct, i.e., from small wave-numbers to 
large ones (see, e.g., Biskamp 2003), and consequently, the 
injection scale can be roughly identified with the system's 
structural scale. In the classical Kolmogorov theory, the 
scale of the largest eddies, which contain most of the en- 
ergy, is referred to as the integral scale. This is also the 
length scale appearing in the definition of the Reyonlds 
number and serves as as a measure for the scale over which 
the turbulent fluctuations are correlated. To avoid compli- 
cations we shall practically identify the integral scale with 
the correlation length and the injection scale, denote them 
as £, and refer to them interchangeably. This identification 
has an obvious intuitive physical basis. 

In a disk that is strongly non-isotropic, two disparate 
injection scales seem a priori to appear - the horizontal 
one -^h^ ''0 and the vertical one ero (both given here 
in dimensional units). It is, however, reasonable to identify 
£ with the smaller of the two, because this is the size of 
largest eddies. A quantitative measure of £ can be "phe- 
nomenologically" determined by the relevant wave-number 
at which the disturbance kinetic energy E]^{k) spectrum 
peaks. The fact that the scale at which resistivity and vis- 
cosity come into play need not be equal (Pm = Rem/Re is 
not generally of order 1) also makes the determination of 
the dissipation scale not unique. In accretion disks, typi- 
cally Pm < 1, i.e., fcj/ > kri where fc^ and krj, respectively, 
represent the dissipation scale wave-numbers for viscosity 
and resistivity. In thin accretion disks, especially if the 
ionization is not complete, these inequalities are strong. Q 
Phenomenologically, the inertial range is considered to be 
the k region, in which E]^{k) exhibits a power-law behavior. 
Thus, the inertial range in the MHD turbulence of a disk 
can be considered as including wave-numbers roughly in 
the range <C ^ fc^. This may be further complicated 
by the fact that MHD turbulence can be inherently non- 
isotropic, even if the system is structurally isotropic (the 
mean magnetic field breaks the symmetry) and one should 
consider fc|| and fc_L separately (this is the case, for exam- 
ple, in the Goldreich & Sridhar theory of interstellar MHD 
turbulence) For a fairly detailed discussion of the MHD tur- 
bulence length scales and relevant references see Biskamp 
(2003). 

If a local (SB) calculation is performed, the scale of the 
system is the box scale, which is (5 <C e for SSB and ^ ~ e 
for LSB. In this context, even though it is evident that the 
equations resulting from the SB approximation possess a 
radial (i.e. in the shear-wise direction x) scale invariance. 



^ Note that in cold disks it is expected that Pm "C 1 through- 
out the bulk of the structure. 



computations done on the equations resulting from the SB 
approximation clearly do introduce a length scale, which is 
the size of the box. Thus, the spatial scale invariance sym- 
metry of the governing equations is lost in the solutions 
whenever a computation is performed. Moreover, any anal- 
ysis done within the SB approximation can be meaningful 
only if the length scale of the physically relevant structure 
(e.g. the most unstable linear mode, the injection scale) is 
significantly smaller than the box size (which we shall re- 
fer to as L, in dimensional units). In particular, nonlinear 
analysis (e.g., performed by numerical simulation) will not 
be able to capture possible nonlinear coupling between the 
unstable linear modes with modes whose wavelength scale 
is larger than L (see also the discussion at the end of Sect. 
[5]). It appears that Goldreich & Lynden-Bell (1965) were 
well aware of this fact and indeed in their nonlinear analy- 
sis they returned to the full equations of motion in rotating 
coordinates. 

Local linear stability analysis (see Balbus & Hawley 
1998) performed on a thin {ho/ro = e <C 1) Keplerian disk 
with disturbances assumed to have the form oc exp(ik • x) 
(i.e., actually having no structure) predicts for the max- 
imally unstable MRI mode • = Uoy/TE/A, with 
'Va = Bo/-\/47rp being the Alfven velocity defined here in 
dimensional units. For a vertical background field we get 
for the vertical wavelength of the most unstable mode, A^, 

27r _ Sir Va _ Sir Va ^ _ Stt ^ _ e Svr ^ 

Here, the thin disk relation ho — Vs/^a a-nd the definition 
of the box size L = S ■ ro have been used. 

In order to capture the fastest growing mode in the 
SSB {S <C e) the plasma beta has to be very large, explic- 
itly /3<; 40e^/(5^. Sure enough, modes with vertical scale 
smaller than Am are also unstable and therefore any nu- 
merical calculation within the SSB scheme with L < Am 
and with periodic-BC (see below) should display the mode 
whose wavelength is equal to the box size as the most un- 
stable one. Indeed, typical SB simulations (see Hawley & 
Balbus 1989) clearly display this feature (as a transient), 
with the CM (also referred to as "channel flows", see also 
below) breaking up later on, nonlinearly. If the SB size, 
or alternatively /?, is not large enough, so as to allow for 
the capture of the entire unstable part of the spectrum, 
it cannot be a priori expected that the simulation yields 
reliable results. In the LSB scheme (with the size of the 
box of the same order of magnitude as the disk thickness), 
the most unstable linear mode will be contained in the box 
for reasonable values of (3. But, a calculation of this sort is 
more involved - vertical structure as well as compressibility 
have to be included. The inclusion of dissipative effects, i.e. 
finite Re, Rcm, further changes the structure of the disper- 
sion relation and a judicious choice of these parameters may 
"push" the unstable part of the spectrum down to within 
the box size, but the large majority of SB calculations for 
the fixed- flux case have been with ideal MHD. 

In any case, possible nonlinear interactions between the 
short vertical wavelength unstable modes, with modes (pos- 
sibly also horizontal ones) whose scale is larger than the SB 
scale, cannot be captured in a SB simulation. It is conceiv- 
able that such interactions may play a role in the instability 
saturation and the development of activity. In the turbulent 
dynamo problem such an interaction between scales is seri- 
ously considered (e.g., Cattaneo & Tobias 2005). We shall 
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also see in Sect. [S] that when one examines a perturbation 
consisting of a wave-packet (rather than just one mode) the 
situation becomes quite complex, even in linear analysis. 

Moving now to a discussion of the small (dissipative) 
length scales, we note that the resistive and viscous scales 
in accretion disks are extremely small in comparison to the 
box scale, even in the SSB case. A numerical SB calculation 
faithfully resolves scales between the system's (box) size 
and some lower limit Znumi the latter of which results from 
the finite accuracy of the numerical scheme. It is easy to 
see that in all existing SB simulations ^num is much larger 
than any realistic dissipation scale (for an accretion disk). 
In ideal MHD calculations, the scale ^num clearly plays the 
role of the dissipation scale, but it is quite doubtful that 
such numerical dissipation faithfully reproduces the correct 
physical dissipation, and so it may conceivably influence the 
results on larger scales. This issue will be discussed in more 
detail in Sect. [H 

We conclude this discussion of relevant scales by point- 
ing out again, and giving some examples of, the above men- 
tioned apparent inconsistency in some SB numerical simu- 
lations that use neither SSB nor LSB, but a kind of "inter- 
mediate" approximation between them. Balbus, Hawley & 
Stone (1996) include compressibility but no vertical struc- 
ture (LSB has to include both, and SSB none), while Sano 
& Stone (2002) include compressibility and investigate fur- 
ther physical ingredients (the Hall effect), but also take the 
background structure as vertically uniform. Similarly, the 
simulations of Papaloizou & Fromang (2007) and Fromang 
et al. (2007) use the LSB (box size equal to the disk thick- 
ness), but the vertical background is uniform. In all these 
calculations periodic-BC are employed in the vertical di- 
rection. Moreover, the initial conditions of the perturba- 
tions in the last two calculations involve a sinusoidal ver- 
tical field in the shear-wise direction x with an amplitude 
that makes the magnetic energy density in the perturbation 
not infinitesimal with respect to the pressure. It means that 
this is a case where finite-size (nonlinear) perturbations are 
causing a subcritical transition to turbulence. This is be- 
cause under these viscous and resistive conditions there is 
no exponentially growing normal mode instability but there 
is, instead, strong transient growth (cf. ROP07). It is not 
clear how simulations considering a medium whose vertical 
structure deviates substantially from that of a disk, have 
insufficient resolution in the ideal case (as shown by the 
authors themselves), and employ periodic-BC in the verti- 
cal, can faithfully uncover the extremely intricate processes 
that must be operating in such transitions (cf. ROP07). We 
think that it would be fair to say that, though the results 
are interesting and worth studying in their own right (at 
least from the standpoint of understanding the interplay of 
the mechanisms at work), the applicability and relevance 
of these results to the physical conditions characterizing 
real disks is doubtful. In contrast, the simulations of Lesur 
& Longaretti (2007) are done in the SSB model (incom- 
pressible flow, constant background) and periodic-BC in all 
three directions. This is consistent according to our view of 
SSB vs. LSB expressed above, but the periodic-BC, when 
applied in a SB of too limited a size, may still pose a diffi- 
culty (see below in Sect. [5]), and it would thus be advisable 
to experiment with different boundary conditions as well. 

It is not yet completely clear how the inconsistencies 
and assumptions discussed above, in the models used thus 
far, may influence the results and how they will ultimately 



relate to real thin accretion disks (or at least large sec- 
tions of them). However, given the concerns we raised thus 
far we think that these matters must be openly reconsid- 
ered. Certainly a lot more comparative numerical experi- 
ments are necessary before a better understanding can be 
achieved. It is our opinion that if the goal is to analyze the 
effects of compressibility and vertical structure in a sec- 
tional representation of a real thin accretion disk then the 
LSB should be used together with a physically meaning- 
ful prescription for the boundary conditions in the vertical 
direction. The LSB by its very derivation includes these ef- 
fects in a physically and mathematically consistent way. In 
contrast, if the SSB (having a uniform background density 
and pressure) is used, there is no need to simulate a com- 
pressible flow, which is computationally more costly and 
physically more involved (e.g., it has to be monitored to 
insure the smallness of the Mach number), even though it 
is not incorrect, per se. 

4. Numerical resolution and its relationship to 
dissipation and saturation 

Astrophysical fluid systems and their magneto-fluid coun- 
terparts (like an accretion disk flow) are endowed with ex- 
tremely large Reynolds (Re) and magnetic Reynolds (Rem) 
numbers and therefore numerical simulations cannot yet 
resolve the full dynamical range of such turbulent flows. 
In other words, given the present magnitude of computing 
power, numerical resolution of the full spatial range of these 
systems, from the system's energy injection scale through 
the full inertial scale and down to the dissipation scale, is 
still out of reach. 

It has been argued in the past (e.g., Balbus & Hawley 
1998) that perhaps a full dynamic range is actually not 
needed in the accretion disk problem because the scale of 
the most unstable linear MRI modes is large and, in the sat- 
urated turbulent state, most of the energy resides in, and 
the angular momentum transport is done by, the large scale 
eddies. Moreover, it has been explicitly stated in that paper 
(basing the argument on the numerical work on hydrody- 
namic turbulence of Oran & Boris 1993) that the physical 
nature of the turbulence on large scales will not change if 
one lets "numerical dissipation" cut off the eddy energy 
spectrum on the grid-scale ^num, instead of allowing for the 
full inertial range, down to the true physical dissipation 
scale. 

This conjecture is debatable for hydrodynamical flows 
(see, e.g., Pouquet, Rosenberg & Clyne 2002), even for Inum 
only modestly larger that the dissipation scale and cer- 
tainly for conditions appropriate for accretion disk flows. 
For an extensive, up-to-date, discussion on the use of unre- 
solved hydrodynamical codes for the study of turbulence, 
see Grinstein, Margolin & Rider (2007). In any case, the 
conjecture appears to be clearly unfounded for MHD tur- 
bulence with very large Reynolds numbers. Boldyrev & 
Cattaneo (2004) found that numerical or experimental in- 
vestigations of MHD turbulence (at least for small P^, as 
is the case in most disks) need to achieve extremely high 
resolution to describe magnetic phenomena adequately and 
Schekochihin et al. (2004, 2005, 2007) have addressed in de- 
tail the question of the relevant scales for the various cases 
(values of Re and Pm) and the resulting computational 
challenge (the need for extremely high resolution). These 
works deal with the turbulent dynamo problem, which is in 
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principle related to the accretion disk problem (at least in 
the zero-flux case) . Similar general conclusions also appear 
in Biskamp (2003). 

Thus, it appears that a viable estimate of a in a MHD 
turbulent disk cannot be obtained from poorly resolved SB 
simulations. In addition, it is even difficult to quantitatively 
assess the energetics of MRI driven turbulence from such 
simulations. Gardiner & Stone (2005) investigated this is- 
sue in detail (sec also Sano & Inutsuka, 2001) using an 
ideal compressible SB calculation with fixed-flux. The dis- 
sipation and hence the thermal energy was calculated by 
relying on "numerical" viscosity and resistivity (Sano & 
Intsuka 2001 used an explicit prescribed resistivity). If the 
nature of the turbulence on scales resolved by these calcula- 
tions (i.e, larger than /num) does indeed depend on the unre- 
solved small scales (see above), the energetics and transport 
properties should be affected as well. The numerical investi- 
gations of Pessah, Chan & Psaltis (2007) and of Fromang & 
Papaloizou (2007) actually show, for the fixed-flux case that 
increasing the resolution of simulations relying on numer- 
ical dissipation shows a decreasing trend in the calculated 
transport, suggesting that this is indeed the case. 

The saturation of the linear MRI (a supercritical transi- 
tion to turbulence in the fixed-flux case), must be effected 
by dissipation if the resulting state of the SB is to rep- 
resent a thin Keplerian disk. Saturation, in general, can 
happen cither by the "removal" of the cause of the linear 
instability and/or by the dissipation of the energy driving 
the linear instability. For the MRI "removal" would be ei- 
ther the expulsion (or redistribution) of the background 
magnetic flux, or the neutralization of the destabilizing 
shear profile. Beginning with the latter we note that the 
gravitational field of the central object, which causes the 
Keplerian shear in the disk, can obviously not be annihi- 
lated by any internal processes in the disk. Neutralization 
of the shear profile can only occur if a new azimuthal ve- 
locity field emerges whose profile masks the destabilizing 
effect of the background shear. The aggregate flow field 
could then be such that the fluid is then marginally stable 
to the MRI. However, this would imply a significant devia- 
tion from a Keplerian profile and the resulting state of the 
system would no longer represent a rotationally supported 
thin disk because there would now emerge significant radial 
gradients of the pressure. With regards to the former satu- 
ration option we note that the surface integrated magnetic 
flux in the SB remains a constant. Therefore, setups consid- 
ering a constant-vertical field cannot vertically expel that 
uniform field (processes resulting from buoyant instabilities 
of the horizontal flelds inside the disk are still possible, as 
well as possible redistribution of flux). Thus, a fully turbu- 
lent state of MRI-driven turbulence must achieve saturation 
by disposing of the energy gained during the linear instabil- 
ity stage by dissipation at the smallest scales, if it is meant 
to represent the saturated turbulent state of a thin disk. 
It is thus difficult to see how an ideal MHD simulation 
that does not resolve the dissipative scales can faithfully 
represent the process of turbulent dissipation and angular 
momentum transport in such a MHD-turbulcnt state of a 
fixed-flux disk. The work of Fromang & Papaloizou (2007) 
suggests that this is the case in a zero-flux system as well. 
Thus proper treatment of the microscopic dissipation seems 
to be imperative. 

Numerical modeling of turbulent flows is obviously not 
exclusive to astrophysics and a large body of work and lit- 



erature exists on this subject (e.g.. Pope 2000, Davidson 
2004, and references therein). Among the basic approaches 
are direct numerical simulation (DNS) and large eddy sim- 
ulation (LES), where in the former approach all scales (up 
to the numerical domain size, of course) arc calculated 
(and well resolved) employing some finite Reynolds num- 
bers, while in the latter approaches scales below Znum are 
treated with the help of a turbulence model. The books 
by Wilcox (1998), Dubois, Jaubereau & Temam (1999), 
Pope (2000) and Davidson (2004) discuss in detail these 
approaches and their subtleties as well as other turbulence 
models. Problems encountered in engineering applications, 
as well as those addressed by the fluid-dynamical commu- 
nity in general, opt for the appropriate numerical scheme, 
according to the problem at hand. The guiding principle in 
making the choice of the method is usually the desire to 
introduce the minimum amount of complexity while cap- 
turing the essence of the relevant physics. In astrophysical 
flows in general and in the accretion disk problem in partic- 
ular, the choice of the BC introduces an additional diSiculty 
for faithful numerical modeling. 

We flnd it useful to classify the following three different 
numerical approaches to the problem of nonlinear develop- 
ment of MHD turbulence in accretion disks. 

1. Local ideal and dissipative calculations - with the pur- 
pose of extracting the maximum information on the lo- 
cal properties of the turbulence in magnetized rotating 
shear flows and obtaining some quantitative measure 
of the turbulent transport coefficients (in particular). 
These are a DNS-like, maximally resolved, calculations. 

2. Global calculations employing some physically moti- 
vated sub-grid scale turbulence model. Such models can 
fully include mass transfer as well as other important 
properties of accretion disks. In principle, these allow 
confronting accretion disk models with relevant obser- 
vations. These are a LES-like calculations. 

3. Global calculations, albeit with unrealistically low 
Reynolds numbers, - in the purpose of discovering the 
trends of the dynamics with increasing Re and/or Re^. 
In a way, this approach is intermediate between the for- 
mer two. 

All existing SB simulations nominally belong to the first 
class, however the great disparity between the smallest re- 
solved scale Inum and the true dissipation scales prevent 
them from being regarded as true DNS. There exist also 
other studies that address the nonlinear development of the 
MRI by utilizing rotating magnetic plane Couette (rmpC), 
whose equations are formally equivalent to the SB equa- 
tions, or magnetic Taylor-Couettc (mTC) setups. These 
are flows conflned between "walls" on which specific well- 
defined boundary conditions (wall-BC) are employed and 
they have been mainly aimed at laboratory experiments of 
the MRI. In some of these calculations spectral schemes 
with very high resolution have been used (e.g., Obabko, 
Cattaneo & Fischer 2006) so as to achieve a true DNS sta- 
tus for relatively very large values of the Reynolds numbers. 

At present, existing global or quasi-global (LSB) calcu- 
lations all suffer from insufficient resolution to achieve the 
status of a DNS simulation of accretion disks with "astro- 
nomical" Reynolds numbers. Except for the old a prescrip- 
tion for MRI mediated or induced MHD turbulence, there 
is as yet no thoroughly tested sub-grid model as might be 
employed in a LES (or one using some other turbulence 
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model) calculation. The MHD turbulence model developed 
by Ogilvie (2003) and the model for hydrodynamic turbu- 
lence in Couette flow (as in the recent work by DubruUe 
at al. 2005 on strongly rotating flows) appear promising. 
Thus far, such models (or any other for that matter) have 
not been employed in simulations pertaining to accretion 
disks. Moreover, turbulence models are usually tailored to 
a particular problem, exploiting relevant results of labora- 
tory and/or numerical experiments. The former model used 
the results of SB simulations, which we discuss in this pa- 
per and doubt their reliability, while the latter model is 
probably inappropriate for MHD flows. 

Among such global calculations Gammie & Balbus 
(1994) investigated the linear stability in an ideal LSB (with 
periodic-BC in the shear-wise direction and various choices 
for the vertical BC) and found that the BC play an impor- 
tant role. Stone et al. (1996) performed a numerical simu- 
lation for these conditions, obtaining a magnetically dom- 
inated disk with a corona and estimated the value of the 
effective a parameter. However, it is difficult to confidently 
evaluate this result in light of our previous comments - in 
particular because the numerical resolution seems insuffi- 
cient and periodic-BC are used in the vertical direction. 
Hawley (2001) performed a global numerical calculation on 
a cylindrical domain, but used a particular artificial viscos- 
ity prescription and it is difficult to say much about the 
effective Reynolds numbers in this calculation (except that 
they are low and uncontrolled). 

The work of Kersale et al. (2004), who investigate a 
global linear problem for a model cylindrical incompress- 
ible flow testing the effects of a variety of BC, including 
those that bring about unstable wall-modes, can be cate- 
gorized as belonging to the above third approach. In the 
follow-up work Kersale et al. (2006) use a spectral code to 
examine the nonlinear development of the wall-modes by 
adopting a well defined (if, however, relatively low) Re at 
P„i = 1. This approach has at least the advantage that 
dissipation is controlled and there may be hope that both 
the linear and nonlinear trends of important flow proper- 
ties can perhaps be asymptotically uncovered for realistic 
Reynolds numbers (like, e.g., in the local hydrodynamic 
study of Lesur & Longaretti, 2005). 

Given the doubt discussed above that numerical dissi- 
pation by itself is sufficient in an LES calculation, it seems 
to be imperative to look for a sub-grid model that is appro- 
priate for the problem at hand. Investigating MRI induced 
or mediated turbulence in local DNS simulations may cer- 
tainly offer in-roads into a better understanding of its local 
salient physical features, which include, e.g., some char- 
acterization of the emerged turbulence in terms of eddy 
spectra and correlations. Information of this sort can pro- 
vide the ingredients for the needed "turbulence model," 
which would go into a faithful global accretion disk LES 
(or a simulation employing other turbulent flow computa- 
tional schemes). Since local simulations are needed to infer 
the turbulent behavior down to the dissipation scale, rmpC 
with physically definite boundary conditions seem to be not 
less appropriate than SB (whose periodic-BC introduce ad- 
ditional difficulties and inconsistencies, see below) for this 
purpose. Of course, such simulations must be fully resolved 
(DNS) and pushed to the maximum Re and Re^ values in 
order for them to yield this kind of information. 

The present limitations on computing power and archi- 
tecture currently rule-out making progress through global 



DNS simulations of accretion disks. Global simulations with 
unrcalistically small Reynolds numbers are a possibility, 
but it is still unlikely that one can learn much about the tur- 
bulent angular momentum transport occurring under real 
accretion disk conditions. However, it is a step closer in that 
direction and it seems (to us) to be a worthwhile tract to 
follow this strategy. 

5. Symmetry and boundary conditions 

The SB approximation with periodic-BC is endowed with 
a particular important symmetry: it allows for shear-wise 
(in the Cartesian geometry of the box it corresponds to the 
x-direction) invariant solutions. The CM are solutions of 
the linear problem in the fixed-flux case and they are un- 
stable (exponentially growing) for wave-numbers kz below 
some critical value for a given = Bq. Moreover, these 
linear modes individually are also exact solutions to the 
nonlinear SB equations, that is, they can appear and ex- 
ponentially grow (on the rotation time-scale) beyond any 
bound. Goodman & Xu (1994) showed that these growing 
modes become linearly unstable and eventually break up. 
This work has served as being an important step toward un- 
derstanding the transition from linear MRI to turbulence in 
accretion disks in the flxed-flux case (see Balbus & Hawley 
1998). 

The i-invariance symmetry, which allows for the exis- 
tence of CM, is obviously broken by any non-periodic-BC 
on X even in local SB analyses and this is manifested on 
the box scale. As discussed in Sect. [3J curvature terms in 
swirling flows (like accretion disks or mTC) as well as ra- 
dial structure gradients, which also break the above sym- 
metry, are global. The SB scale, appropriate for thin accre- 
tion disks satisfies S ■ vq < L < e ■ tq, while radial symmetry 
breaking by the above mentioned global physical effects are 
manifested on a length scale of the order tq. This symmetry 
breaking can thus be "communicated" through a real com- 
pressible disk flow only on timescales larger than the hor- 
izontal sound-crossing time, Tsound ~ rg/vs ~ (ro//io)/f^o- 
On the other hand the (linear) growth time of the CM is 
TcM ^ 1/f^o- Thus, because ro//io = ~ 50— 100 in most 
thin accretion disks, one might think that the broken sym- 
metry can not affect the development of the CM. However, 
the saturated MHD turbulent state in the disk, postulated 
to result from the Goodman & Xu (1994) CM instability, 
must persist for very long times (relative to even Tgound) 
and so this global state has to be driven by a continuous 
nonlinear process, giving rise to energy injection along the 
full global scale of the disk. 

It is important to mention, in this context, the work 
of Coppi & Keyes (2003), who showed that in a disk of fi- 
nite vertical extent, in which the modes have to localized 
vertically (and not an infinite cylinder or a finite one, but 
with periodic BC), the modes necessary acquire a charac- 
teristic oscillatory profiles in the radial direction, i.e., are 
certainly not of the CM type and have much longer rise- 
times. This is probably related to the very recent investi- 
gation by Liverts & Mond (2007), who examined whether 
in the simplest ideal MHD setting, with the only spatial 
dependence of the perturbations being oc e''^^^ (like of the 
CM), an absolute or perhaps convective linear instability 
(see, e.g., Schmid & Henningson 2000, §7.2, for a detailed 
discussion of this topic) occurs. Solving analytically an ini- 
tial value (Cauchy) problem for the MRI in the fixed-flux 
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case, with an initially localized perturbation in the form 
of a Gaussian with prescribed width, they found that the 
spatio-temporal evolution of the packet is not as trivial as 
that of single CM or a superposition thereof. The initial per- 
turbation was found to split into two wave packets moving 
in opposite directions (vertically) from the initial position. 
As each wave-packet propagates at its group velocity (of 
the order of the Alfven speed), the amplitude of the pack- 
ets does not grow in time. Instead, it oscillates sinusoidally 
in time in the co- moving frame of the wave-packet. In addi- 
tion to the two oppositely moving wave-packets, Liverts & 
Mond found an absolutely unstable singular component, re- 
sulting from the well-known branch point of the dispersion 
relation (a;f,,fcf,) = (3i/4, \/l5/4), growing in time asymp- 
totically (that is for large enough t) like cx e^* j \ft^ where 
7 = \ijJh\ and having an amplitude that is a function of 
the initial packet width. Consequently, the linear stability 
analysis yields a situation considerably more complicated 
than just of exponential growth of CMs, with their subse- 
quent secondary instability of Goodman & Xu (1994). In 
particular, since the disk is not infinite in its vertical ex- 
tent, reflections of the wave packets by the boundaries may 
complicate the situation. The absolutely unstable compo- 
nent does not grow fast enough before the reflected packets 
return to their origin and therefore periodic-BC are simply 
inapplicable. Thus, the quite simple and elegant Goodman 
& Xu local linear scenario, which is generally accepted as 
an explanation of the mechanism for transition from linear 
MRI to turbulence in a fixed-flux disk, should be reconsid- 
ered. 

Linear analysis of a global configuration does not show 
the presence of CM (see the discussion of "body modes" in 
Kersale et al. 2004). Moreover, the nonlinear development 
also appears to be quite different (see Kersale et al. 2006). 
Only a global calculation, of the kind performed by Kersale 
et al. (2006), for example, but allowing for a compressible 
flow and definite vertical BC, can perhaps clarify this point. 
In any case, it seems that the role of the CM in transition 
to global MHD turbulence in an accretion disk threaded by 
a vertical magnetic field is quite academic. As important 
as such local understanding may be, it is not clear if it has 
any relevance to the global problem. 

Turning now to a discussion of the effects of periodic- 
BC, we first note again that the majority of numerical stud- 
ies on the nonlinear development of MHD turbulence in ac- 
cretion disks, aiming at the understanding of the properties 
of the saturation and transport, employed some version of 
the SB approximation. The periodic-BC employed in vir- 
tually all SB calculations, have several adverse side-effects. 
We shall distinguish here between the periodic-BC on x in 
sheared coordinates and the vertical (on z) periodic-BC. 

It was already mentioned in the previous section, that 
the z periodic-BC may cause the inability to capture the 
most unstable mode. But this can be remedied by either 
making /3 very large (i.e. a very small magnetic field) or by 
including dissipation and choosing an appropriate combi- 
nation of Re,Pm, and /SQ (see Umurhan, Menou & Regev 
2007). However, z periodic-BC for the fixed-flux case are 
bound to cause un-physical results, for at least two reasons. 
First, they would give rise to a spurious reentry of any con- 
vected component of the perturbation, when it reaches a 

^ Or, alternatively, the the box Cowling number, Co = 
Vi/illL^ = 1/13 



boundary (Liverts & Mond 2007, see above). Second, mag- 
netic fields, as is well-known, have a necessarily global na- 
ture. This means that any local calculation, in which well- 
ordered magnetic field lines extend beyond the computa- 
tional domain, must employ BC based on the knowledge 
of the magnetic configuration outside the domain. Thus, 
arbitrary vertical BC (in particular z periodic-BC) on the 
field (in the fixed-flux case) will necessarily exert extrane- 
ous magnetic stresses on the system and may substantially 
alter its dynamics. This point has recently been brought to 
attention by Shu et al. (2007). 

We finally turn to the most important observation about 
the use of periodic-BC in numerical studies of turbulence 
and its bearing on SB simulations. Clearly, space periodic- 
BC are not accessible in realistic physical situations, and 
they only can be faithfully employed in the study of homo- 
geneous turbulence when one assumes that real boundary 
effects are not important (Dubois, Jauberteau & Temam, 
1999). The caveats pertaining to the use of periodic-BC 
in numerical simulations of turbulence are discussed in de- 
tail in §7.2 of Davidson (2004) from which we only bring 
here issues that are explicitly relevant to SB simulations 
of MHD turbulence in an accretion disk. In order for the 
bulk of the turbulence not to be seriously affected by the 
imposed periodicity, the box size should satisfy L ^ £, oth- 
erwise the "tails" of the relevant correlation functions (of 
the turbulent fluctuations) are not only cut off, but also 
very significantly lifted (see figure 7.6 of Davidson 2004). 
This crucially affects the behavior of the large scale dynam- 
ics. Thus, SSB is clearly out of the question if the aim is to 
gain any understanding of accretion disk turbulence up to 
the injection scale. If we accept the reasonable conjecture 
that the injection scale relevant for an accretion disk is in- 
deed £ = ivf^ ho, the minimal size of the SB must be the 
disk height. But then z periodic-BC cannot be used and 
we are only left with the option of semi-global (LSB) or 
global simulations. SSB simulations can still be useful, but 
only for understanding the dynamics on very small scales, 
with the purpose of devising an appropriate sub-grid model 
for global, or at least semi-global LES. In the next section 
we shall demonstrate, with the help of a numerical experi- 
ments, that using a SB, whose size is not large enough (as 
compared to the scale of the relevant dynamical structures), 
and periodic-BC, gives rise to spurious energy fluctuations, 
casting a serious doubt on the physical significance of the 
simulation results. 



6. Energetics of the SB 

6.1. Integral statements 

In any continuum mechanics problem posed in a finite do- 
main, as is the case for SB MHD numerical calculations 
of the kind discussed here, it is generally advantageous to 
consider integral physical conservation statements, valid on 
the domain. A particularly useful consideration of this kind, 
in this problem, arises from the examination of the energy 
budget for the system of equations ((T][5]). Using the defini- 
tions B = {bx,by, Bo + bz}, the total magnetic field vector 
and U = {u,Uq + v,w}, the total velocity vector, where 
the background held Bq may be equal to zero (for the zero- 
flux case) and the the background shear is embodied in 
Uq = —qriox, one may take the appropriate inner products 
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with the equations of motion and integrate them over the 
SB. 

After rather standard manipulations the following total 
energy theorem results 



in which 



(10) 



(11) 



is the integral over the SB of the total energy density 
£(x, t), defined to be 



hear 



-EB-q%x', (12) 



where Eu = ^jup is the kinetic energy density of the per- 



turbation atop the linear shear Uo, Eb = jpB the total 
magnetic energy density, i^shcar = —qfloxv + ig^figX^ and 
we have explicitly used Uq = —qflox. The term —qfl^x^ in 
the definition of £ is the centrifugal potential energy den- 
sity, arising in a SB rotating with the local Keplerian veloc- 
ity (referred to as "tidal" -potential by Gardiner & Stone, 
2005). It obviously merely represents a local Taylor expan- 
sion of the gravitational field potential upon which the fun- 
damental shear in this fiow derives its existence. 

The terms on the righthand-side of equation (fTU]) are 

S EE - ^ n • IJ{£ + vj)d(7 + ^<j>h- B(B • U)dcr -t- 

-f^/(A.V)U^.. + ^/(n.V)B^d., (13) 

where fi is the unit normal of the bounding surface of the 
box and da is the differential area element of the bounding 
surface, and 



1 

Re 



1 



Re. 



(14) 



where we have used the notation |VUp = X^i l^^iP for 
i = x,y, z and similarly for B. Clearly, the surface term 
S (jl3p includes possible energy injection into the domain 
by influx of matter as well as through viscous and resistive 
stresses. The term T) represents bulk dissipation. 

6.2. Numerical experiments 

We have performed a number of numerical experiments 
with the purpose of examining the effects of the kind of 
BC employed (periodic versus wall) and of the SB size (rel- 
atively to a relevant integral scale £ of the problem) on the 
evolution of the SB energy content. We have chosen to fo- 
cus on the total energy £ of the domain, as this quantity 
has (at least to us) the most definite and best understood 
physical meaning. 

The essence of the problem, as we see it, is embodied in 
the BC applied in the shear-wise (x) direction and there- 
fore we opt for simulating the evolution of the SB dynami- 
cal equations, restricted to two dimensions - x and y. The 
vertical dimension, as important as it may be (we have dis- 
cussed, for example, the problems that arise from periodic- 
BC in z) should have no bearing on the effects we wish to 



demonstrate here. Using a 2-D code makes the calculation 
tractable, but, as is well-known, one cannot expect a per- 
sistent activity in two dimensions as any 2-D turbulence 
ultimately decays if it is not maintained by some kind of 
forcing. In a 3-D SB simulation the activity need not a priori 
decay, but this fact should not be important for our basic 
findings, which comprise of effects occurring on a short (as 
compared to the viscous decay) time scale. We shall fur- 
ther comment on this matter below, after describing our 
calculation and its results. 

We shall consider the evolution of the total domain en- 
ergy for two sets of boundary conditions. The first of these 
are the periodic-BC, as laid out before in the paper. The 
second type of BC are appropriate for a rotating channel 
and we refer to them as wall-BC They implement at the 
shear- wise boundaries i) no normal flow; ii) no normal mag- 
netic fleld flux; and iii) free-slip boundary conditions; i.e., 



u — 0, bx = 0, dxVy — 0, at x = ±Lx/2. 



(15) 



Wall-BC are periodic in the azimuthal direction. 

Note that the equations solved are identical for both 
sets of BC - they are the two-dimensional restriction of the 
SSB equations ([l])-(l5]). We take (3—1 and remark that 
since there is no z dependence in our 2-D setup, Bq (the 
background vertical field) does not play any dynamical role 
in the problem. In addition, it is convenient in this case to 
use the flux function $ and the stream function f/;, such 
that 



dxi) = V, dyip = -u, 



and 



9.* 



by, dy^^-bx, (16) 



and so the solenoidal property of the velocity and magnetic 
field disturbances is automatically insured. It is perhaps 
instructive to also remark that the energy integral (I10|) ap- 
propriate for our 2-D numerical calculation takes on the 
form 



e^qno 



uv - 



bxby 



Lx/2 



dy,-V. 



for periodic-BC, and 



Re 



v{--q^lo) 



-ix/2 



dy-V 



(17) 



(18) 



for wall-BC. Thus, as is apparent, the surface term S re- 
duces to a simple and readily understandable form. 

The code we use in the experiment with periodic-BC 
is a modification (so as to include the MHD terms) of the 
spectral evolver we developed for the purely hydrodynamic 
problem, see Umurhan & Regev (2004) for details. The nu- 
merical experiment employing wall-BC also consists of a 
spectral code, but there a different spectral function basis 
is used, appropriate for wall-BC. In both cases, we evolve 
the flux function <I> and the stream function tp, instead of 
u, V, bx and by, as explained above. We temporally evolve 
the resulting set of equations using the modifled Crank- 
Nicholson scheme as implemented in Umurhan & Regev 
(2004). We investigate the dynamics of a flow with Re = 
2000 and Rm = 700 and make sure that the resolved dy- 
namical scales are in the viscous and resistive regimes so 
as to guarantee that the dissipation is fully resolved in all 
simulations. For both types of BC we run the simulations 



for a SB of size Lx = n and Ly 
of double size (in Lx , that is 



2tt as well as for an SB 
27r). 
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The resolution in the x-direction in the small domain 
(La; = TT and Ly = 2tt) is taken to be = 64 and Uy = 
64. Because of our experience with externally driven shear 
problems (e.g. Umurhan & Regev, 2004), we are careful 
to to increase the resolution in the shear-wise direction on 
account of the strong crenellation that the external shear 
creates out of the disturbances (the Orr-mechanism) , so 
that for the double (in the shear- wise extent of the domain) 
case, we double the resolution too, i.e., when = 2tt we 
take = 128. 

The simulations are initiated with white noise in the 
vorticity (lu — V^i/') the source J of the flux function 
(J = V^$). The initial uj and $ are shown in Fig. [TJ This 
initial disturbance is identical in all simulations and it is 
localized away from the radial boundaries. We do this to 
gain better control of the effects that the boundaries have 
in the evolution. The only relevant spatial dynamical scale 
that may be perceived as the integral, or injection, scale is 
the extent of the initial perturbation, suggesting that £ 2, 
initially. 

In the upper panel of Fig. [2l we show the evolution 
of the total energy in the SB of shear-wise extent — 
TT, for periodic-BC and wall-BC, while this figure's lower 
panel depicts the same for a SB of double shear-wise extent 
(Lx = 2tt). The fluctuations in £ depicted on the upper 
panel in periodic-BC simulations are quite dramatic and 
we have found that they correlate with the action taking 
place on the boundaries. Clearly, in the smaller SB we have 
(, that is we are in the situation in which the periodic 
box size is too small, as we have discussed at the end of 
the previous section. The dynamics, for the periodic-BC 
case, is seriously afi^ected by the artificial distortion of the 
correlation function brought about by the periodic-BC on 
a too small domain. In contrast, the evolution with wall- 
BC proceeds smoothly, as the boundary conditions do not 
allow artificial energy inflow/outflow. The fluctuations in 
the periodic-BC case have a time scale related to the period 
time of the shearing box given by 

T = 



Fig. 2. The total energy (in arbitrary units) in the SB as 
a function of time (in units of ^Iq ). In the upper panel, 
the size of the SB is La; = tt and Ly — 27r and the result of 
the calculation with periodic-BC is compared with the one 
using wall-BC. In the lower panel, the same is shown for 
a SB having a double shear- wise extent. In both cases the 
simulations are started with the same initial conditions, as 
depicted in Fig. [1] and the evolution is shown during the 
period in which the typical (of 2-D) decay of the activity 
has not yet fully set in. 



as we can clearly see from the time behavior of the surface 
term (not shown). This behavior is also apparent in Fig. 
[21 but is less prominent and clear. We reason that these 
fluctuations are driven by successive passages of the de- 
veloped coherent structures (vortices, in this case) in the 
imaginary neighboring boxes. The vortices that appear in 
all our simulations are actually somewhat tighter than the 
initial vorticity perturbation extent, suggesting that £ ~ 1 
or so. We also see this trend played out in the lower panel, 
however an important difference is clearly apparent. The 
domain size is now 6, that is, considerably larger than 
£. In this case, the periodic-BC do not induce very signifi- 
cant spurious energy exchange with the exterior and conse- 
quently the fluctuations of the SB energy content are rather 
mild. We may reasonably conjecture that a calculation with 
Lx ^ £ would reveal no difference between the pBC and 
wBC curves. 

The longtime decay, which we have mentioned above as 
being typical in a 2-D problem of this kind, is obviously 
due to the bulk dissipation embodied in the term T). We 
comment again that the effects we wished to demonstrate 
in the above numerical experiment can not be significantly 
altered by the third (vertical) direction because the key 
factor responsible for them is the boundary condition in 
the shear-wise direction. 

In summary, our numerical experiments show that the 
type of boundary conditions, coupled with the box size, 
adopted in investigations of local disk turbulence have im- 
portant implications on the physical results. When one con- 
siders the total energy of disturbances within the formalism 
of the SB equations, we note that the total energy 2; flue- 
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tuates dramatically on account of the injection/extraction 
of energy across the boundaries when the size of the box 
is not large enough as compared to the relevant injection 
scale. It is thus impossible to sec how a SB simulation with 
periodic-BC can reveal any dynamical behavior that is rel- 
evant to a true accretion disk, in which I has to be at least 
of the disk thickness size. 

7. Discussion and Summary 

Based on the arguments of the previous sections we are 
led to conclude that existing SB calculations, done with 
periodic-BC in the fixed-flux as well as zero-flux cases, suf- 
fer from some serious inconsistencies and problems. These 
inconsistencies make these calculations quite unreliable for 
their purpose - finding the mechanism of and quantitatively 
estimating the turbulent angular momentum transport and 
energy dissipation in accretion disks. As we have discussed 
in this paper, the problems stem from the limitations of the 
box scale, insufficient resolution of the numerics, unjustified 
symmetry properties (due to periodic-BC) and, as an ad- 
ditional consequence of periodic-BC and limited box size, 
the possibility of unrealistic energy sources (sinks), which 
can seriously interfere with the dynamics. 

In the fixed-flux case the linear MRI gives rise to a su- 
percritical transition. In this case the SB calculations with 
periodic-BC manifest channc;l flows, resulting from the non- 
linear exponential growth of the CM, whose problematic 
nature arises not only from their shear-wise symmetry but 
also the fact that the linear instability appears to be more 
involved than just the one comprising of a linear super- 
position of the CM. The zero-flux simulations in the same 
SB setup do not allow for CM, but they suffer from all 
the other difiiculties. Indeed, also zero-flux numerical cal- 
culations, when tested for convergence (most recently those 
of Fromang & Papaloizou 2007) , indicate that in the ideal 
case a (a measure for the angular momentum transport) 
decreases steadily with increasing resolution. Moreover, in 
the zero-flux case a dynamo action must be invoked and the 
MRI is thought to play a role in it. But it is not clear how 
such an intricate self-sustained process, in which a subcriti- 
cal transition, similar to the one occurring in hydrodynamic 
shear flows (see ROP07), can be faithfully uncovered by SB 
periodic-BC simulations, with their problems, as discussed 
in this paper. 

KPL07 conjectured that the fixed-flux case is proba- 
bly irrelevant to accretion disks, but clearly accretion disks 
exist that are threaded by external fields (like, e.g., those 
around young stellar object, see Shu et al. 2007). Moreover, 
Pessah, Chan & Psaltis (2007) claim the opposite, i.e. that 
MRI related angular momentum transport can only occur 
in disks threaded by strong (~ equipartition value) fields, 
basing this conclusion on their numerical simulations. In 
any case, even if the ultimate development of the mag- 
netic field in the disk makes its geometry inside very in- 
volved, this fact cannot annihilate the external field whose 
lines must close through the disk. It seems, therefore, more 
likely that the apparent almost linear dependence of a on 
Bz in some SB fixed-flux simulations is simply due to the 
non-reliability of these calculations. Moreover, SB zero-flux 
calculations also suffer from most of the same difficulties as 
the fixed- flux ones (as discussed above). 

The most recent numerical studies using the SB with 
periodic-BC in both fixed-flux and zero-flux case have in- 



corporated explicit Re and Re™ (that is, physical non- 
zero viscosity and resistivity) and studied the dependence 
of the angular momentum transport (or a) on these non- 
dimensional numbers. As we have indicated before, this ap- 
proach can yield valuable information. Lesur & Longaretti 
(2007) and Fromang et al. (2007) have reported that a oc 
PJ5j with k; > 0. We have shown earlier, using asymptotic 
semi-analytical techniques, (Umurhan, Menou & Regev, 
2007) that such scaling can be naturally expected for P^ <C 
1, and also found that k depends on the boundary condi- 
tions (URM07). Our work was done for flxed-flux conditions 
and in a thin-gap mTC setup (which eliminates the shear- 
wise invariant CM). However, in general, any dependence 
of the angular momentum transport on a positive power 
of Pto makes the role that the MRI takes as the driver of 
angular momentum transport in astrophysical disks with 
such astronomical Reynolds numbers rather doubtful (at 
least in the fixed-flux case, where it is considered to be the 
sole driver). 

Laboratory experiments have so far not yielded any con- 
clusive indications of MRI driven turbulence, save for the 
case in which relatively strong azimuthal fields are applied 
(HoUerbach & Riidiger 2005; Stefani et al. 2006; Szklarski 
& Riidiger 2007). In thin accretion disks, however, strong 
toroidal (or any other, for that matter) magnetic fields 
cause significant magnetic pressure and are likely to be 
inappropriate for accretion disks in the thin, rotationally 
supported, disk limit. 

It is not easy to answer the question, what kind of nu- 
merical scheme or which setup is best suited for the prob- 
lem? In any case, it is not clear what is less appropriate 
for local numerical studies of accretion disks: wall-BC or 
periodic-BC in the SB. We believe that the use of periodic- 
BC, especially in the radial direction, suffers from serious 
problems unless the power in the perturbation fields is lo- 
calized away from the boundaries of the system (Dubois et 
al. 1999, Davidson 2004). 

We think that if one wants to extract reliable phys- 
ical information on the local nature of MHD turbulence 
in which shear is instrumental, on scales from the dissipa- 
tion one and up to well into the inertial range, it is likely 
that mrpC (e.g., ROP07) or mTC (e.g., Liu, Goodman & 
Ji 2006; Obabko, Cattaneo & Fischer 2006) in a thin-gap 
setup are more appropriate. There certainly are no walls 
in accretion disks, but wall-BC calculations of this kind 
can, at least, be made mathematically consistent and phys- 
ically sound. One may think that periodic-BC (which are 
also not realistic for accretion disks) are better suited, but 
as we have shown SB calculations with periodic-BC suffer 
from difficulties and inconsistencies for this problem, math- 
ematical as well as physical. 

It seems to us that the logical approach to the accretion 
disk angular momentum transport problem should thus 
consist of: 

1. mathematically and physically sound calculations, most 
likely in the mrpC or mTC thin-gap setups, for the 

purpose of learning about the detailed local properties 
of MHD turbulence (its energy spectrum, correlations 
etc.); 

2. devising a turbulence model, based on the above, so as 
to replace the a model; 

3. global calculations employing the above turbulence 
model(s) and various boundary conditions, giving rise 
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to accretion disk models that can be confronted with 
observations. 

Regarding the first item, it has to be kept in mind that the 
saturation mechanism may perhaps be different from the 
case when there are no physical boundaries (see Knobloch 
& Julien 2005, who have performed an asymptotic MRI 
analysis, but for a developed state, far from marginality) . 
The hope is, however, that this should not affect too much 
the local properties of the MHD turbulence. 

The previously mentioned approach of performing 
global calculations with low Reynolds numbers (e.g., 
Kersale et al. 2005) are also useful as complementary stud- 
ies, in which the dependence of the results on relevant non- 
dimensional numbers can be uncovered. In any case, any nu- 
merical study must be controlled, in the sense that one can 
faithfully (that is, in a numerically resolved way and with- 
out introducing any a-priori periodicity of disturbance cor- 
relations) experiment with the relevant mathematical prop- 
erties (e.g., BC, symmetries) and physical quantities (e.g., 
viscosity, resistivity) and uncover their effects. Schemes in- 
cluding hyper-viscosity, which are routinely used in numer- 
ical turbulence studies, may be employed, but this non- 
trivial tool must be used with proper care. Global quanti- 
ties (like the energy) must obviously also be monitored and 
their development fully understood. 

The problem of angular momentum transport in accre- 
tion disks, especially if MHD turbulence in a non-isotropic 
strongly sheared medium is the physical mechanism respon- 
sible for it, is extremely difficult and involved. This work, 
as well as that of KPL07 and others, indicates that the un- 
derstanding gained so far from SB and other numerical sim- 
ulations seems unfortunately not to be adequate enough to 
provide a viable physical prescription that can be effectively 
used in modeling accretion disks, beyond the old a- viscosity 
parametrization. It seems that only an extensive and coher- 
ent program, of the kind that has been undertaken in the 
turbulent dynamo problem (see Schekochihin et al. 2007), 
is likely to lead to significant progress. In such a program 
a massive numerical effort is indispensable, but we would 
like to stress, in concluding, that analytical (asymptotic or 
other) approaches can also be very useful in guiding the 
numerical and laboratory work and better understanding 
the underlying physics. 
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